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A novel simulation model of cloud microphysics is developed, which is named Super-Droplet 
Method (SDM). SDM enables accurate calculation of cloud microphysics with reasonable cost 
in computation. A simple SDM for warm rain, which incorporates sedimentation, condensa- 
tion/evaporation, stochastic coalescence, is developed. The methodology to couple SDM and a 
non-hydrostatic model is also developed. It is confirmed that the result of our Monte Carlo scheme 
for the coalescence of super-droplets agrees fairly well with the solution of stochastic coalescence 
equation. A preliminary simulation of a shallow maritime cumulus formation initiated by a warm 
bubble is presented to demonstrate the practicality of SDM. Further discussions are devoted for the 
extension and the computational efficiency of SDM to incorporate various properties of clouds, such 
as, several types of ice crystals, several sorts of soluble/insoluble CCNs, their chemical reactions, 
electrification, and the breakup of droplets. It is suggested that the computational cost of SDM 
becomes lower than spectral (bin) method when the number of attributes d becomes larger than 
some critical value, which may be 2 ~ 4. 



I. INTRODUCTION 

Although clouds play a crucial role in atmospheric phe- 
nomena, the numerical modeling of cloud is not well es- 
tablished up to the present. The macro-scale processes 
such as the fluid motion of moist air associated with 
clouds is called "cloud dynamics" and the micro-scale 
processes such as condensation/evaporation, stochastic 
coalescence, and sedimentation of water droplets, are 
called " cloud microphysics" . These two processes mu- 
tually affect each other, and thus we can say that 
cloud formation and precipitation development are typi- 
cal multiscale-multiphysics phenomena. Cloud dynamics 
model to describe the fluid motion in the atmosphere 
has been well developed (Jacobson 2005). However, it is 
still difficult to perform an accurate simulation of cloud 
microphysics though several simulation methods, such 
as bulk parameterization method (Kessler 1969; Ziegler 
1985; Murakami 1990; Ferrier 1994; Meyers et al. 1997), 
spectral (bin) method (Berry 1967; Soong 1974; Bott 
1998, 2000), and the exact Monte Carlo method (Gille- 
spie 1975; Seefielberg et al. 1996), have been proposed. 
Numerical methods to accurately simulate both processes 
and their interactions are required to understand and pre- 
dict cloud-related phenomena. 

We have developed a novel, particle-based simula- 
tion model of cloud microphysics, named Super-Droplet 
Method (SDM), which enables accurate numerical simu- 
lation of cloud microphysics with reasonable cost in com- 
putation. Though several extensions and validations are 
still necessary, we expect that SDM provides us a new ap- 
proach to the cloud- related open problems, such as the 
cloud and aerosol interactions, the cloud-related radia- 
tive processes, and the mechanism of thunderstorms and 
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lightning. 

SDM can be regarded as a sort of Direct Simula- 
tion Monte Carlo (DSMC) method, which was initially 
proposed to simulate the Boltzmann equation for pre- 
dicting rarefied gas flows (Bird 1994). Particularly, 
SDM has similarity to the Extended version of the No- 
Time-Counter (ENTC) method, which was developed 
by Schmidt and Rutland (2000) for the simulation of 
spray flows. Though the basic idea of using a compu- 
tational particle with varying multiplicity is common in 
both SDM and ENTC, there is a difference in the proce- 
dure of the Monte Carlo scheme. 

The main purpose of the present paper is to de- 
velop a simple SDM for warm rain, which incorporates 
sedimentation, condensation/evaporation, and stochas- 
tic coalescence, and discuss the theoretical foundations 
of SDM. The methodology to couple SDM and the non- 
hydrostatic model is also developed. Further discussions 
on the extensions and the computational efficiency of 
SDM are also carried out, even though briefly. 

The organization of the present paper is as the follow- 
ing. In Sec. II we introduce the principle model, which is a 
very primitive microphysics-macrodynmics coupled cloud 
model and a good starting point for our discussions. In 
Sec. Ill we review the traditional methods, i.e., the ex- 
act Monte Carlo method, bulk parameterization method, 
and spectral (bin) method and see how they are derived 
from the principle model and what are the problems they 
involve. In Sec. IV, we develop a simple SDM for warm 
rain, which can be regarded as a coarse-grained model 
of the principle model, and discuss the theoretical foun- 
dation and characteristics of SDM. Numerical simulation 
scheme is also proposed for each process. Especially, a 
novel Monte Carlo scheme for the stochastic coalescence 
process is developed and validated. In Sec. V a very pre- 
liminary result of a shallow maritime cumulus formation 
initiated by a warm bubble is presented to demonstrate 
the practicality of SDM. In Sec. VI we briefly discuss the 
possibility to extend SDM to incorporate various proper- 
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ties of clouds, such as, several types of ice crystals, several 
sorts of soluble/insoluble CCNs, their chemical reactions, 
electrification, and the breakup of droplets. The esti- 
mation of the computational cost of SDM suggests that 
SDM will be computationally more efficient compared to 
spectral (bin) method when the number of attributes d 
becomes larger than some critical value, which may be 
2^4. Finally, a summary and concluding remarks are 
presented in Sec. VII. 



II. THE PRINCIPLE MODEL 

This section is devoted for the introduction of the prin- 
ciple model. The principle model is a very primitive 
microphysics-macrodynmics coupled cloud model. The 
characteristic point of this model is the very precise de- 
scriptions of cloud microphyscs processes, which is re- 
alized by following all the water droplets in the atmo- 
sphere. Though being not so useful because of the high 
cost in computation, this model gives us a clear overview 
of the relationship between each methods to be intro- 
duced later, i.e., the exact Monte Carlo method, bulk 
parameterization method, spectral (bin) method, and 
super-droplet method. 



A. cloud microphysics for the principle model 

1. definition of real- droplets 

Warm clouds and rain consist of large number of water 
droplets with broad range of radius from micrometers to 
several millimeters. Several sorts of aerosol particles are 
also floating in the atmosphere, which can be a seed of a 
cloud droplet. These particles are called Cloud Conden- 
sation Nuclei (CCN). We regard each of these particles in- 
cluding both water droplets and CCNs as a real-droplet. 

Let N r (t) be the number of real-droplets floating in 
the atmosphere at time t. Each real-droplet is located 
at position Xi(t), i — 1, 2, . . . , N r (t), with velocity Vi{t) 
and radius Ri(t). Soluble CCN is dissolved in each real- 
droplets. Here, we consider only one soluble substance 
as the CCN and denotes the mass of solute dissolved 
in the real-droplet i. We refer to these state variables, 
the velocity Vi(t), the radius Ri(t), and the mass of so- 
lute Mi(t), as the attributes of real-droplets denoted by 
at (t) . Note that we do not regard the position Xi (t) as an 
attribute for convenience. Assuming that the droplet ve- 
locity immediately reaches their terminal velocity, Vi(t) 
is determined by their radius Ri(t). As a result, the num- 
ber of independent attributes of our real-droplet is two, 
the radius Ri (t) and the mass of solute Mi (t) . In the fol- 
lowing, we give the time evolution law of real-droplets. 



2. motion of real-droplets 

In general, the motion equation of a real-droplet is de- 
termined by the drag force from the ambient air and the 
gravitational force: 



dvj 
1 dt 



dx ■ 

mig + F D (vi,U(xi),Ri), — p = v u (1) 

dt 



where rrii := (An / 3) py lq is the mass of the droplet i 
with pii q = 1.0 g cm~ 3 , g is the gravitational constant, 
Fd is the atmospheric drag, and U(xi) is the wind ve- 
locity at the droplet position Xi(t). Assuming that all 
the droplets reach their terminal velocity immediately 
enough, Vi(i) is evaluated as Vi(t) = U(xi) — zVooffii), 
which is the steady solution of fT]). We give v 0o (Ri) 
according to the semi-empirical formulas developed by 
Beard (1976). Consequently, Vi(t) is no longer an inde- 
pendent attribute of real-droplet in our model. 



3. condensation/ evaporation of real- droplets 

We consider the mass change of real-droplets through 
the condensation/evaporation process according to 
Kohclcr's theory, which takes into account the solution 
and curvature effects on the droplet 's equilibrium vapor 
pressure (Koheler 1936; Rogers and Yau 1989). The 
growth equation of the radius Ri is derived as the fol- 
lowing. 
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Here, S is the ambient saturation ratio; Fk represents the 
thermodynamic term associated with heat conduction; 
Fd is the term associated with vapor diffusion; the term 
a/ Ri represents the curvature effect which expresses the 
increase in saturation ratio over a droplet as compared 
to a plane surface; the term b/Rf shows the reduction 
in vapor pressure due to the presence of a dissolved sub- 
stance and b depends on the mass of solute Mi dissolved 
in the droplet. Numerically, a ~ 3.3 x 1CP 5 /T (cm) and 
b ~ A.3iMi/m s (cm 3 ), where T (K) is the temperature, 
i ~ 2 is the degree of ionic dissociation, m s is the molecu- 
lar weight of the solute. R v is the individual gas constant 
for water vapor, K is the coefficient of thermal conductiv- 
ity of air, D is the molecular diffusion coefficient, L is the 
latent heat of vaporization, and e s (T) is the saturation 
vapor pressure. 



4- stochastic coalescence of real-droplets 

Two real-droplets may collide and coalesce into one 
big real-droplet and this process is responsible for pre- 
cipitation development. It is worth noting here that the 
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coalescence process is of our particular concern in this pa- 
per because SDM is expected to overcome the difficulty 
in the numerical simulation of coalescence process. 

Our primary postulate, which is commonly used in 
many previous studies, is that the coalescence process 
can be described in a probabilistic way. Consider a re- 
gion with volume AV which is small enough compared 
to the space-scale of the atmospheric fluid. Then, we 
assume that the real-droplets inside this region are well- 
mixed and coalesce with each other with some probability 
during a sufficiently short time interval (t, t + At), i.e., 
there exists such C(Rj,Rk) that satisfies 



Pjk =C(Rj,R k )\vj - v k 



At 
AV 



'-K{Rj,Rk) 



At 
AV 



(3) 



^probability that real-droplet j and k 
inside the small region AV will coalesce 
in the short time interval (t, t + At). 

Here, C(Rj, R k ) can be regarded as the effective collision 
cross-section of the real-droplets j and k and may be 
evaluated as 



C(Rj,R k ) = E(Rj, Rk)ir(Rj + Rk) 2 , 



(4) 



where E(Rj,Rk) is the collection efficiency, which takes 
into account the effect that a smaller droplet would be 
swept aside by the stream flow around the larger droplet 
or bounce on the surface (Davis 1972; Hall f980; Jonas 
1972; Rogers and Yau 1989) . K(Rj,R k ) is defined by 
K(Rj,Rk) := C(Rj, Rk)\vj —Vk\ and it is called the coa- 
lescence kernel, which will later appear in the stochastic 
coalescence equation. 

The cloud microphysical process of the principle model 
has been defined. We will turn to determine the cloud 
dynamics process of the principle model. 



Here, D/Dt := d/dt + U ■ V is the material derivative; 
p = p c i + p v is the density of moist air, which is repre- 
sented by the sum of dry air density p e i and vapor den- 
sity p v ; q v — Pv/ P is the mixing ratio of vapor; U is the 
wind velocity; T is the temperature; 9 is the potential 
temperature; n = (P/ P )( R <i/ c p) is the Exner function 
with reference pressure Pq; p w is the density of liquid 
water; S v is the source term of water vapor associated 
with condensation/evaporation process; g is the gravita- 
tional constant; A and k are the transport coefficient; R^ 
is the gas constant for dry air; c p is the specific heat of 
dry air at constant pressure; L is the latent heat of vapor- 
ization. Equations ([5])- ([9]) represent equation of motion, 
equation of state, thermodynamic equation, mass conti- 
nuity, and water continuity, respectively. Note here that 
it is assumed that q v <C f to derive these equations. 

There are three coupling terms from the microphysics: 
p w represents the momentum coupling; LS v /c p H is 
the release of latent heat through the condensa- 
tion/evaporation process; S v is the source term of vapor 
through the condensation/evaporation process. These 
source terms are evaluated by the microphysics variables, 
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where rrii is the mass of real-droplet i. 

The principle model has defined completely now and 
we can see that the description of the cloud microphysics 
is very precise. In the next section, we review the tra- 
ditional methods from the point of view how they are 
related to the principle model. 



III. TRADITIONAL METHODS 



B. cloud dynamics of the principle model 

We may apply non-hydrostatic model to describe the 
cloud dynamics process of the principle model. There are 
several versions of non-hydrostatic model, but the basic 
equations used in this paper is as the following. 

P^r = -VP-(p + Pw )g + XV 2 U 1 (5) 
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In this section, we review the traditional methods and 
see how they are derived from the principle model and 
what are the problems involved. 



A. The exact Monte Carlo method 

The exact Monte Carlo method was developed by 
Gillespie (1975) and improved by SeeBelberg et al. 
(1996) . This method simulates the principle model very 
directly. Their procedure repeatedly draws a random 
waiting time for which the next one pair of real-droplets 
will coalesce. Consequently, it is very rigorous, but not 
so useful to simulate a cloud formation process in a suf- 
ficiently large region because of the very high cost in 
computation. 
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B. bulk parameterization method 



T2l) reduces to 



Bulk parameterization method is most commonly used 
to deal with cloud microphysics. This method empirically 
evaluates the source terms from microphysics, such as p w 
and S v , as a function of the field variables in the cloud 
dynamics model (Kessler 1969; Ziegler 1985; Murakami 
1990; Ferrier 1994; Meyers et al. 1997). Then, the gov- 
erning equations become a closed form in cloud dynamics 
variables and we need not cope directly with cloud micro- 
physics. Consequently, the number of degrees of freedom 
is reduced so much that the computational demand is ex- 
tremely relaxed. However, the accuracy seems to be not 
always enough. This is very likely because bulk param- 
eterization is justified under the assumption that there 
exists a closed- form equation in cloud dynamics, though 
this is not guaranteed. 



C. spectral (bin) method 

Another approach to cloud microphysics is spectral 
(bin) method (Berry 1967; Soong 1974; Bott 1998, 2000). 
Spectral (bin) method is characterized by solving the 
time evolution equation of the number density distribu- 
tion of real-droplets. 

Let Xi — (47r/3)i?f be the volume of the real-droplet i, 
and n(X, M, x,t) be the number density distribution of 
real-droplets, i.e., n(X,M,x,t)AXAMAV denotes the 
number of real-droplets with their volume in the range 
of (X, X + AX), mass of solute in the range of (M, M + 
AM), inside the space region AV around x, at time t. 
Then, the time evolution equation for n(X, M, x, t) which 
corresponds to the principle model can be derived as the 
following. 
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dt ~ 2 
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where v t (X,x) := U(x) — zv^R) is the terminal ve- 
locity of a droplet, f g (X,M,x) := AttR x (r.h.s of ©) 
satisfies dX/dt = /„, K is the coalescence kernel defined 
in ©, X" := X - X', and M" := M - M' . The l.h.s 
of ([P2]) corresponds to the advection of real-droplets in 
real-space and attribute-space, and the r.h.s corresponds 
to the expected dynamics of the stochastic coalescence 
process. 

On the ideal case when there is only the stochastic 
coalescence process and neither sedimentation nor con- 
densation/evaporation process exists, and further, only 
the droplet volume X is regarded as the attribute, then 



dX'n{X')n(X")K(X',X") 



n(X) / dX'n(X')K(X,X'). 
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This is the simplest case of spectral (bin) method and 
([T3")) is the equation which is sometimes called the 
stochastic coalescence equation (SCE), but in this paper, 
SCE is used as a general term refers to the governing 
equation of spectral (bin) method. 

These equations are derived from the probability 
law © under the assumption that the number density 
has no correlation (Gillespie 1972). This assumption 
seems to be correct (Seefielberg et al. 1996), thus the 
solution of SCE can be very accurate. However, the com- 
putational cost to solve SCE is very expensive, especially 
if we consider many types of microphysics processes and 
the number of attributes d becomes large. We will come 
back to this remark in Sec. VI, but in a word, this is 
owing to the fact that the governing equation contains a 
d-multiple integral. 



IV. SUPER-DROPLET METHOD 

In this section, we develop a simple SDM for warm 
rain, which can be regarded as a coarse-grained model 
of the principle model, and discuss the theoretical foun- 
dation and characteristics of SDM. Numerical simulation 
scheme is also proposed for each process. Especially, a 
novel Monte Carlo scheme for the stochastic coalescence 
process is developed and validated. 



A. cloud microphysics for Super-Droplet Method 

1. definition of super-droplets 

First of all, let us define what a super-droplet is. 
Each super-droplet represents a multiple number of real- 
droplets with the same attributes and position, and the 
multiplicity is denoted by the positive integer £i(t) € 
{1,2,...}, which value can be different in each super- 
droplet and time-dependent due to the definition of co- 
alescence introduced later. That is to say, each super- 
droplet has their own position Xi(t) and own attributes 
at (t) , which characterize the £j (t) number of same real- 
droplets represented by the super-droplet i. Here, the 
attribute di(t) consists of the radius Ri(t) and the so- 
lute mass Mi(t) for this simplest case. Note here that 
no two real-droplets have exactly the same position and 
attributes. In this sense, super-droplet is a kind of 
coarse-grained view of real-droplets both in real-space 
and attribute-space. Let N s (t) be the number of super- 
droplets floating in the atmosphere at time t. Then, these 
super-droplets represent N r (t) = X^i^i City) number of 
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real-droplets in total. In the following, we give the time 
evolution law of the super-droplets. 



2. motions of super- droplets 

Except the coalescence process, each super-droplet be- 
haves just like a real-droplet. Thus, a super-droplet 
moves according to the motion equation ([5]) and Vi(t) 
is evaluated by the terminal velocity. 

It is easy to simulate this process numerically. Define 
the time step A.t rn and evaluate the terminal velocity 
Vi(t) and the position X{(t + At m ) = Xi(t) + At m Vi(t) at 
each time step. 



3. condensation/ evaporation of super- droplets 

Condensation/evaporation process of super-droplets is 
also governed by the growth equation . 

We may adopt implicit Euler discretization scheme for 
the numerical simulation of yielding 



before coalescence 



E>2 _ p2 

2At„ 



(5-1) 



Ri, 



+ 



w. 



F k + F d 



where Ri, n ~ Ri(t — nAt g ). It may be not possible to de- 
termine analytically the next step value Ri >n+ \ from the 
previous value Ri. n . However, using Newton- Raphson 
scheme, can be evaluated numerically. 



4- stochastic coalescence of super- droplets 

Stochastic coalescence process for the super-droplets 
must be formulated carefully so that the super-droplets 
well-approximate the behavior of real-droplets. The for- 
mulation procedure is in two steps: 1) define how a pair 
of super-droplets coalesce; 2) determine the probability 
that the super-droplet coalescence occurs. 

Consider the case when a pair of super-droplets (j, k) 
will coalesce. These super-droplets represent £j number 
of real-droplets (a,j,Xj) and number of real-droplets 
(a,k,Xk). Let us define that exactly min(£j,£fc) pairs 
of real-droplets will contribute to the coalescence of the 
super-droplet pair (j,k). Figure. [T] is a schematic view 
of the coalescence of super-droplets. In this example, 
£j = 3 and = 2. We can see that min (£,j,£,k) = 2 pairs 
of real-droplets undergo coalescence, which results in the 
decrease of multiplicity ^ : 3 — > 1 and the increase of the 
size of the super-droplet k. 

Below is the complete definition of how a pair of super- 
droplets (j, k) change their state after the coalesce: 

1. if £j 7^ £fc, we can choose £j > without losing 



^=3 



c 



after coalescence 



£=2 




1=1 



FIG. 1: Schematic view of the coalescence of super-droplets. 
Two super-droplets with multiplicity 2 and 3 undergo coa- 
lescence (upper left). This represents the coalescence of 2 
real-droplet pairs (lower left and right). As a result the super- 
droplet with multiplicity 2 becomes larger and the multiplicity 
of the other super-droplet decreases 3 — > 1 (upper right). 



generality, and 

£j = £j _ £fc j Ck = 6s 
= Rj, R' k = (Rj 
M'j=Mj, M'j — (Mj + Mfc), 



*2) 1/3 , 



x ■ 
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(14) 
(15) 
(16) 
(17) 



where the dashed valuables represent the updated 
value after the coalesce. 
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(18) 

(19) 
(20) 
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where Gauss' symbol [x] is the greatest integer that 
is less than or equal to x. 

Note here that in the case £j = no real-droplets will 
be left after the coalescence because all £j number of real- 
droplets contribute to the coalescence. Thus, we divide 
the resulting £j number of real-droplets into two super- 
droplets. 

This definition of super-droplet coalescence possess has 
a favorable property that the number of super-droplets 
is unchanged in most cases though the number of real- 
droplets always decreases. The number of super-droplets 
is decreased through the coalescence only when £j = = 
1, i.e., both super-droplets are real-droplets. This results 
in £'j = and = 1, and we remove the super-droplet j 
out of the system. Because the number of super-droplets 
corresponds to the accuracy of SDM, the number conser- 
vation of super-droplets suggests the flexible response of 
SDM to the drastic change of the number of real-droplets. 
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Having defined how a pair of super-droplets coalesces, 
we turn to determine the probability that the coalescence 
occurs. If we require the consistency of the expectation 
value, the probability is determined as the following. 



=probability that super-droplet j and k 
inside the small region AV will coalesce 
in the short time interval (t, t + At c ). 



(22) 



Indeed, we can confirm that the expectation value is 
exactly the same to that of the principle model (real- 
droplet view). The super-droplet j represents £j num- 
ber of real-droplets with attribute and super-droplet 
k represents number of real-droplets with attribute 
cifc. From the real-droplet view, there are number 
of real-droplet pairs which have the possibility to coa- 
lesce with the probability Pjk, i.e., the number of real- 
droplet pairs which will coalescence follows the binomial 
distribution with £j£k trials and success probability Pjk- 
Thus, from the real-droplet view, the expectation value 
of the number of coalesced pairs is £,j£,kPjk- From the 
super-droplet view, a coalescence of super-droplets j and 
k represents the coalescence of min(£j, £&) pairs of real- 
droplets with attribute aj and a^. Thus, the coales- 
cence of super-droplets j and k is expected to represent 
min (Ci>'?fc)-Pjfc ) = min(^,4)max(^,^)P j7£ = / ',/. 
number of coalescence of real-droplet pairs, which is equal 
to the value of the real-droplet view. It is worth notic- 
ing here that the variance of super-droplet view becomes 
larger than that of real-droplet view. From real-droplet 
view, the variance of the number of coalesced pairs is 
^jtkPjkO- - Pjk) — £j£kPjk =■ V r (Poisson distribu- 
tion limit). From super-droplet view, the variance is 

{wmfe,£ k )} 2 P$ - {tjtkPjk} 2 - min(0,6)K, which 
is min(£j, times larger than that of the real-droplet 
view. 

We have defined the cloud microphysics of super- 
droplets. If all the multiplicity = 1 then our model is 
equivalent to the principle model. Thus, we can change 
the accuracy of our model by changing the initial dis- 
tribution of multiplicity = 0)}, or equivalently, the 
initial number of super-droplets N s (t — 0). If the con- 
vergence to the principle model is rapid, we can say that 
the use of SDM is beneficial. 

The definition of the coalescence process of super- 
droplets is rather theoretical and we have to do numerical 
simulation on this model somehow. In the next section, 
we develop a Monte Carlo scheme for the stochastic co- 
alescence, which can simulate the coalescence process of 
super-droplets efficiently. 



5. Monte Carlo Scheme for super-droplet coalescence 

Let us introduce a grid which covers the real-space. 
We can choose any kind of grids, which may be not uni- 
form, but the volume of each cell must be small enough 



compared to the space-scale of the atmospheric fluid so 
that the coalescence occurs according to the probability 
(f2"2"| in each cell. Let us make a list of super-droplets in 
a certain cell at time t: I := t2, ■ ■ ■ « n ,}, where n s is 
the total number of super-droplets in the cell. Then, 
by definition, each pair of super-droplets (j, fc) G I 2 , 

j =/= fc, coalesces with the probability P$ within the 

short time interval (t,t + At c ). Note that P.^ <C 1 be- 
cause At c is small. We have to examine all the possible 
pairs (j,k) £ I 2 , j ^ fc, to know the random realization 
at the time t + At c . However, this yields computational 
cost 0(n 2 ) and this is not efficient. Instead, we propose 
a novel Monte Carlo scheme which leads us to 0(n s ) cost 
in computation. 

Let L := {(J x , ki), (j 2 , k 2 ), ■ ■ ■ , (j[„»/2], *K/2])} be a 
list of randomly generated, non-overlapping, [n s /2] num- 
ber of pairs. Non-overlapping means that no super- 
droplet belongs to more than one pair. This list can 
be made by generating a random permutation of I and 
make pairs from the front, which cost 0(n s ) in computa- 
tion (see the appendix). By examining only these [n s /2] 
pairs instead of the whole possible n s (n s — l)/2 pairs, 
the cost will be 0(n s ). In compensation for this simpli- 
fication, all the coalescence probability is divided by the 
decreasing ratio of pair number, and the corrected, scaled 
up probability for the ath pair is obtained as 



P 



is) 



n s (n s 



jc,k a 



(23) 



This may be justified if the following relation holds good, 



^/2] 



3,k,j^k 



j-k,j^k 

which represent the consistency of the expectation num- 
ber of coalesced real-droplet pairs in this cell. 

Based on the above idea, below is the complete pro- 
cedure of our Monte Carlo scheme to obtain one of the 
stochastic realizations in a certain cell at time t + At c 
from the sate of super-droplets at time t. 

1. Make the super-droplet list at time t in this certain 
cell: I = {ii,«2, ■ ■ - in s }- 

2. Make the [n s /2] number of candidate-pairs 
L = {(ji,ki),(j2,k 2 ),---,(j[ ns /2],k[ ns / 2 ])} from 
the random permutation of /. 

3. For each pair of super-droplets (j Q , k a ) 6 L, gener- 
ate a uniform random number <f) a 6 (0,1). Then, 
evaluate 



7a := 




if 4> a <p a - \p a ] 
if <j) a >p a - \p a ] 



7 



4. If 7 a = 0, the crth pair (j a , k a ) is updated to f+Af c 
without changing their state. 

5. If 7 Q 0, choose £j a > £k a without losing general- 
ity and evaluate 7 Q := min(7 Q , [£j a /£k a ])- 

( a ) ^ _ 7a£fc Q > 0, 
Mj a = M Ja , Mj a = (% M ja + M fc J , 



( b ) If - 7a6o = 0, i.e., 7q 



£• =K*„/2], & = - [&„/2], 



R 



3 \V3 



(%M ja +M ka ), 



X ja ~ X 3" > 



x. 



If £j a = 0, the super-droplet j a is removed out 
of the system. 

Our Monte Carlo scheme numerically solves the 
stochastic coalescence process of super-droplets defined 
in Sec. IV A 4. We used a sort of random sampling of 
candidate-pairs to achieve the computational cost 0{n s ). 
Further extension are implemented in our Monte Carlo 
scheme to deal with multiple coalescence. The integer 7 a 
represents how many times the pair (j Q , k a ) will coalesce, 
and 7q is the restricted value of 7 Q by its maximum value 
[£j a /£k a } ■ Primarily, p a must be smaller than 1 because 
p a represents a probability, thus 7 Q should be either 
or 1. However, if we take At c too much larger, p a may 
become larger than 1. Though the situation 7 Q > 1 is 
not consistent with the fundamental premise of our the- 
ory, we formulate our Monte Carlo scheme to cope with 
this situation so that SDM robustly well-approximate the 
principle model irrespective of the choice of At c . 

Alternatively be said that the most rigorous selection 
criteria of the simulation time step for coalescence process 
At c is that p a <C 1 holds good for all a. Let us estimate 
At c very roughly. The typical number density and radius 
of small cloud droplets may be 10 9 m -3 and R = 10 /Ltm. 
The typical terminal velocity of a rain droplet may be 
1 m s _1 . Then, we may roughly evaluate p a as follows. 



Eir(R ja + R k J 2 \v Ja - v ka \At c 



p a ~ n s P}J ka 

_ n s max , £, ka ) 

^7 ""("j.TLtj iv Ja 

~ 10 9 rrT 3 • 1 • tt(2 x l(T 5 m) 2 • 1 m s _1 • At c 
= 0.47rAt c < 1. 

Thus, it is estimated that At c < 1/OAtt s ~ 0.8 s. It is 
worth noticing here, that our estimation is independent 
of the number of super-droplets n s and the volume of the 
coalescence cell AV. 



6. validation of our Monte Carlo scheme 

In this section, we validate our Monte Carlo scheme 
for the stochastic coalescence process by confirming that 
the numerical results agree well with the solution of the 
SCE ip]). 

The SCE (|13p is the time evolution equation of n(X, t) 
in the ideal system in which only the stochastic coales- 
cence process exists. Correspondingly, the super-droplets 
in this section undergo only the coalescence process ac- 
cording to our Monte Carlo scheme. The super-droplets 
have velocities, but do not move out of the coalescence 
cell, nor do they undergo condensation/evaporation. The 
radius Ri is the only attribute. 

To compare the simulation result, we have to recon- 
struct the number density distribution n(X, t) from the 
data of super-droplets {(£,, Ri)}. There are several ways 
to do this, but we adopt kernel density estimate method 
with Gaussian kernel (Terrell and Scott 1992). Some 
brief discussions on the application of kernel density es- 
timate method to SDM is also developed in Sec. IV B. 
It is convenient to plot the results by the mass density 
function over lni?, which is defined by g(hiR)dhiR := 
pii q Xn(X)dX . The corresponding estimator function 
g(hxR) becomes, 



1 



g{\nR) := —^^miWailnR-lnRi), 



2na 



exp (-Y 2 /2a 2 



Here, AV is the volume of the coalescence cell and 

— 1 / 5 

(j = (TqNs with some constant ag. Theoretically 
the most efficient choice of CTo is estimated as cr — 
(2y/nAV 2 J g" 2 dlni?/M tot 7V r )" 1/5 : wh ere M to t is the to- 
tal mass of liquid water. Because g is the estimated func- 
tion itself, we have to estimate also do from the data of 
super-droplets to find the maximum likelihood estimator 
function g, but in this paper we just choose do empiri- 
cally. 

We have to give an explicit form of the coalescence 
kernel K(Rj,Rk) in ([3]) to determine the dynamics. For 
Golovin's kernel K{Rj,R k ) = b(Xj + Xk), we have the 
analytic solution of SCE (fT3"|) (Golovin 1963) though this 
choice of coalescence kernel is not realistic. The result of 
our comparison simulation for Golovin's kernel is shown 
in Fig. [2k.. We set b = 1.5 x 10 3 s^ 1 . We need only 
one big coalescence cell, the volume of which is chosen 
to be AV = 10 6 m 3 . The time step is fixed as At c — 
1.0 s. The initial number density of real-droplets is set 
to be no — 2 23 m~ 3 . The initial size distribution of the 
real-droplets follows an exponential distribution of the 
droplet volume Xi which is determined by the probability 
density p(X*) = (1/X ) exp (-X t /X ), X = (4^/3)^, 
i?o = 30.531 x 10 -6 m. This setting results in the total 
amount of liquid water 1.0 g m~ 3 . The initial number of 
super-droplet N s is changed as the simulation parameter. 
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The initial multiplicity is determined by the equality 
£i = noAV/N s . cr is fixed as cr = 0.62. 

Figures [3t> and 2c show the results of the case of so- 
called hydrodynamic kernel, which is a much realistic co- 
alescence kernel. Hydrodynamic kernel is defined by ([3]) 
and (J3J. The same coalescence efficiency E(Rj,Rk) is 
adopted as described in Seefielberg et al. (1996) and 
Bott (1998). For small droplets the dataset of Davis 
(1972) and Jonas (1972) is used, and for large droplets 
the dataset of Hall (1980) is used. Because the analytic 
solution is not available for hydrodynamic kernel, we gen- 
erated the reference solution numerically using Exponen- 
tial Flux Method (EFM) developed by Bott (1998, 2000). 
A logarithmically equidistant radius grid is used in EFM. 
We adopted 1000 bins ranging from 0.62 (im to 6.34 cm. 
The time step was set to 0.1 s. These parameter settings 
provide a sufficiently accurate numerical solution of the 
SCE (fl~3|) . In Fig. [2)3, simulation settings are same as case 
(a) except the coalescence kernel. In Fig. f2t we changed 
the initial size distribution by reducing Rq to one-third, 
i.e., i?o = 10.177 x 10~ 6 m, and increased the initial num- 
ber density of real-droplets as n = 3 3 • 2 23 m -3 . Con- 
sequently, the total amount of liquid water 1.0 g m -3 is 
unchanged from case (b). At c = 0.1 s and <tq = 1.5 for 
case (c). 

In case (a) and (b), we can see that the result of SDM 
with N s — 2 13 agrees fairly well with the solution of 
SCE (TIB")) , and much improved when N s — 2 17 . How- 
ever, case (c) is not so good as the previous two cases. 
Even N s = 2 17 is not good and we need N s = 2 21 
for good agreement. This situation typically reveals the 
weak point of SDM. The initial size of the droplets are 
comparatively small and the coalescence hardly occurs. 
We can see that the thick solid line at t = 1200 s is not 
so changed from the initial distribution. However, once 
a few big droplets are created, then the coalescence pro- 
cess is abruptly accelerated. Accordingly, for an accurate 
prediction, we have to resolve the right tail of the distri- 
bution at the time t = 1200 s. However, existence prob- 
ability of super-droplets is very low for such region and 
sampling error occurs if the number of super-droplets N s 
is not sufficient. As a result, we need much more super- 
droplets for case (c) compared to that of (a) and (b). For 
the practical applications to cloud formation, this fact 
would not impose severe restrictions on SDM because the 
condensation/evaporation process dominates the growth 
of small droplets in such cases like (c). 

Anyway, we have confirmed that SDM reproduces the 
solution of the SCE (TIB"]) if the number of super-droplets 
N s is sufficiently large. These results support the validity 
of our Monte Carlo scheme for the stochastic coalescence. 



B. cloud dynamics for Super-Droplet Method 

For the cloud dynamics of SDM, we can also apply the 
non- hydrostatic model (O-©- Note here that we have 
to evaluate the source terms from microphysics (| 1 1) - (|TT|) 
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FIG. 2: Time evolution of the mass density distribution 
g(\nR, t), which is numerically obtained by the Monte Carlo 
scheme of SDM. (a): The case of Golovin's kernel. The solid 
line represents the analytic solution of SCE (|13fl . (b) and (c): 
The case of hydrodynamic kernel. The solid line represents 
the approximate solution of SCE (|13|l . which is numerically 
obtained by EFM. In each case, we can see that the result of 
SDM agrees fairly well with solution of SCE (|13p when the 
number of super-droplets N 3 is sufficiently large. At c = 1.0 s 
for (a) and (b), and At c = 0.1 s for (c). 
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not by real-droplets but by super-droplets: 



>0M) = ^2£,zmi(t)S 3 (x - Xi(t)), 



S v (x, t) 



-1 dpu 
p(x,t) dt 



(24) 



(25) 



For the numerical simulation of (O-Q, we may adapt 
4th-order Runge-Kutta scheme for the time derivative, 
second-order central difference scheme for the spatial 
derivatives, with artificial viscosity to stabilize the nu- 
merical instability. The source terms (fM]) and ([23)1 must 
be evaluated on the numerical grid point of fluid simula- 
tion x*. Generally, introducing an appropriately chosen 
coarsening function w(x), it can be evaluated as 

p w (x*,t) = J w{x* - x)p w (x, t)d 3 x 

i=l 

Here w satisfies J w(x)d 3 x — 1, and would be a piecewise 
linear, localized function with length-scale comparable to 
the grid width. If a super-droplet is simply shared by 
the 8 adjacent grid points at the corner of the cell, then 
w(x) = 1/8 AV inside the cube with volume 8AV, and 
w(x) = outside the cube. 

SDM has been defined completely and the numerical 
simulation schemes have been also proposed. This is the 
main result of the present paper. In the next section, we 
show a preliminary result of SDM. 



V. DEMONSTRATION 

A very preliminary result of a shallow maritime cumu- 
lus formation initiated by a warm bubble is presented in 
this section to demonstrate the practicality of SDM. 

The simulation domain is 2-dimensional {x-z), 12.8 km 
in horizontal direction and 5.12 km in vertical direction. 
Initially, a humid but not saturated atmosphere is strat- 
ified, which is almost absolutely unstable under 2 km in 
altitude. There is a temperature inversion layer above 
2 km. The equations below determines the initial base 
sounding of the atmosphere. 

U{x,t= 0) = 0, 
T(x,z= 0,t = 0) = 305.0 K, 
P(x, z = 0,t = 0) = 101325 Pa, 

<9T f-9.5 K km" 1 if z < 2.0 km, 
~dz~ ~ I +3.0 K km" 1 if z > 2.0 km, 



q v (x,t = 0) = 0.022 exp 



(z + 200.0 m) 
2000.0 m 



The upper and lower boundary is fixed to their initial 
value and the horizontal boundary is periodic. A warm 



bubble is inserted at the center of the horizontal axis 
according to the equation below. 



(x,t = 0) = 9 b {x,t = 0) + 1.0 K 



x exp 



(z - 500.0 m) 
400.0 m 



(a: -6.4 km) V 
1200.0 m J 



where 6b(x,t = 0) is the base state of potential temper- 
ature. 

We simulate the cloud dynamics according to the non- 
hydrostatic model introduced in Sec. II B with the nu- 
merical schemes proposed in Sec. IV B. The grid width 
is Ace = 3.125 m and Az = 4.0 m with time step 
At = 0.0125 s. We set the artificial viscosity v = K = 
1.56 m 2 s _1 . 

Initially, NaCl CCNs are uniformly distributed in space 
with number density 1.0 x 10 7 m- 3 . The solute mass 
Mi(t = 0) has an exponential distribution given by the 
probability density 



/ , ,r \ 1 ( Mi 

p(Mi ) = — exp(- — 



Mn 



1.0 x 10" 16 



Because it is unsaturated at the beginning, the growth 
equation |2]) has a stable and steady solution, which is 
adopted as the initial value of the radius Ri (t = 0). The 
coalescence efficiency used in Sec. IV A 6 is also adopted 
to this simulation. 

For the numerical simulation of the microphysics pro- 
cesses, we also use the numerical scheme proposed in 
Sec. IV. Initially, the super-droplets are uniformly dis- 
tributed in space with number density 1.92 m- 3 with 
common multiplicity, i.e., £i(t = 0) = INT(1.0 x 
10 7 /1.92) = 5208333. The grid for our Monte Carlo 
scheme of stochastic coalescence is same to that of non- 
hydrostatic model. The time steps are At m = At g = 
1.0 s and At c = 10.0 s. 

Cloud started to form at about 8 min and began to 
rain at about 20 min. After raining for half an hour, 
the cloud remained for a while but finally disappeared 
at about 90 min. The amount of precipitation is about 
1 mm in total. Figure [3] is a snapshot at the time 1608 s. 
We plot all the super-droplets with the color and the 
alpha transparency which are determined by their radius 
Ri and multiplicity £j. The color map is indicated in 
the figure and the alpha value is determined by on = 
7.8 x 10- 5 (i?i/10- 3 m) 2 which is proportional to i? 2 ^ 
We can see that there is a turbulent like structures inside 
the cloud. Note that there are super-droplets also in the 
black region, but it is not depicted in this figure because 
the radius is very small, which represent CCNs. 

In this section, we showed a very preliminary result of a 
shallow maritime cumulus formation initiated by a warm 
bubble, and demonstrated that SDM actually works as 
a cloud resolving model. Application to more realistic 
situations are our future works. 
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VI. DISCUSSIONS 
A. extension of Super-Droplet Method 

In Sec. IV, SDM was developed only for warm rain 
with only one soluble substance as the CCN. However, it 
is very conceivable that we can extend SDM to incorpo- 
rate various properties of clouds, such as, several types 
of ice crystals, several sorts of soluble/insoluble CCNs, 
their chemical reactions, electrification, and the breakup 
of droplets. Here, we give a rough overview of how to 
achieve these extensions. 

For this purpose, let us extend the type of attributes 
a. To treat several kinds of soluble/insoluble chemicals 
for the CCNs, new attributes are necessary to represent 
the mass of each chemicals. To treat several types of 
ice crystals, one attribute may represent the state of the 
super-droplets. Additional attributes may represent the 
electric charge and the temperature of the droplet. If 
we do not give a terminal velocity, then the velocity of a 
droplet is also an independent attribute. 

Let us denote the attributes of our new super-droplet 
as a(t) = (a^\ a^ 2 \ . . . , a^), where d is the number of 
independent attributes. Then, we have to give a time 
evolution law for each attributes, which determines the 
individual dynamics of super-droplets: 

^ = f(oi,A(xi)), i = 1,2,..., N s (t), (26) 
at 

where A(xj) represents the field variable which charac- 
terize the state of ambient atmosphere. Equation (|26p 
may include the growth equation ^ and the time evo- 
lution of solute mass dMi/dt = 0. For example, the 
chemical reactions and the electrification process may be 
also incorporated in (f2"6"|) . 

The coalescence process must be specified for our new 
super-droplets. At first, we have to determine what kind 
of real-droplet a' will be created after the coalescence of 
real-droplets dj and This can be either determinis- 
tic or stochastic. Then, the coalescence of super-droplets 
will be determined similarly to (|T¥|) - (f^Tj) . It is not the 
case in (ITTl) and (12"T]) , but we may also change the position 
of super-droplets after the coalescence if necessary. Sim- 
ilar to ([3]), the coalescence probability for real-droplets 
must be given in the form 

P jk =C(aj,a fc )^|«j - v k \ 

^probability that real-droplet j and k 
inside the small region AV will coalesce 
in the short time interval (t, t + At c ), 

which in general is a function of the attributes Oj and 
cifc. Then the coalescence probability for super-droplets 
will be determined by ([22]) . 

The coalescence process is a short time and short range 
interaction between the droplets, but we may also con- 
sider other types of interactions if any. 



Breakup of droplets can be also taken into account. 
This can be done if a governing law is given for the 
breakup process. However, note here that breakup re- 
sults in the increase of super-droplets, and this is not 
preferable because the more the super-droplet is the more 
the computational cost will be. Breakup process may 
result in the explosion of the number of super-droplets 
N s (t) — > oo. In other words, SDM is formulated suitable 
for the simulation of such phenomena in which the coa- 
lescence process is predominant and the breakup process 
has a secondary importance. This is exactly the case for 
cloud formation process. 

SDM can be regarded as a sort of Direct Simula- 
tion Monte Carlo (DSMC) method, which was initially 
proposed to simulate the Boltzmann equation for pre- 
dicting rarefied gas flows (Bird 1994). Particularly, 
SDM has similarity to the Extended version of the No- 
Time-Counter (ENTC) method, which was developed 
by Schmidt and Rutland (2000) for the simulation of 
spray flows. Though the basic idea of using a compu- 
tational particle with varying multiplicity is common in 
both SDM and ENTC, there is a difference in the proce- 
dure of the Monte Carlo scheme. SDM is using randomly 
generated, non-overlapping, [n s /2] number of candidate 
pairs, and allows multiple coalescence for each pair. On 
the other hand, ENTC method chooses the candidate 
pairs randomly from the set of all the possible pairs. The 
number of candidate-pairs to be sampled varies propor- 
tionally to the maximum coalescence probability in the 
cell. Both SDM and ENTC method result in the com- 
putational cost 0(N S ). Probably, ENTC method is also 
applicable to cloud simulations and SDM is applicable to 
spray simulations, but further verifications are still nec- 
essary to compare their capability. 

B. computational cost and accuracy 

Our estimation suggest that the computational cost of 
SDM becomes much lower than spectral (bin) method 
when the number of attributes d becomes large. We 
would like to show the basis of this assertion very briefly. 

Let us estimate the computational cost and accuracy 
of spectral (bin) method. We have extended SDM in the 
previous section, and the corresponding, general form of 
SCE can be derived as follows. 

^M+V..W + V .{/n}=(^) el (27) 

where n(a, x, t) is the number density distribution and / 
is defined in (|26|) . The term (dn/dt) c corresponds to the 
stochastic coalescence process and this will be d-multiple 
integral in general because we have to survey all over the 
rf-dimensional attribute-space to evaluate the variation 
of droplet number density at a. Obviously, this term is 
the most expensive for computation, so we neglect the 
advection terms in the l.h.s of (|27p and focus our atten- 
tion on the simplified equation dn/dt = (dn/dt) c , which 
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corresponds to the ideal case that only the stochastic co- 
alescence exists as the microphysics process. 
We choose Integrated Squared Error (ISE), 



C = d d a {n(a)-n b (a)Y 



(28) 



to evaluate the accuracy of spectral (bin) method. Here, 
rib (a) is the approximate solution generated by spectral 
(bin) method, and C evaluate the difference of n&(a) from 
the exact solution n(a). Let Nb be the number of bins 
for each attribute, i.e., the number of grid points for the 
discretization of n(a) per attribute. Then, if the accuracy 
of spectral bin method is fcth order in attribute-space, 



C scales like C 



N, 



-2k 



Thus, we can estimate the 



amount of operation count and memory needed for the 
computation of spectral (bin) method as the following, 



operation ~ N% d ~ ( — 



memory ~ N d 



2d/k 



d/k 



(29) 



Let us estimate the computational cost and accuracy of 
SDM. To compare the result with spectral (bin) method, 
we only consider the stochastic coalescence process, and 
let us estimate the number density distribution n(a) from 
the super-droplets {(£i,a,i);i — 1, 2, . . . , N s } using the 
kernel density estimation method, which was originally 
developed to estimate the generating probability distri- 
bution from its random sample (Terrell and Scott 1992). 
The density estimator function n(a) with Gaussian ker- 
nel Wa^ (a) is defined by 



i(a) :=J2tM d) (a- ai ), 



(30) 



WW(a) 



-exp{-a 2 /2cr 2 }. 



The error function corresponding to ISE (f2"8")) is the Mean 
Integrated Squared Error (MISE) defined by 



C(a) = E 



d d a {n(a)-h(a)Y 



(31) 



Remember that each {(£t,£ii)} is a random realization 
and C(a) is defined as an ensemble averaged value. Sev- 
eral calculations and assumptions are necessary to derive 
the results, but we can estimate the expectation value 
and variance of n(a) as 



E[n(a)]-n(a)^?- J ]T 



d 2 n{a) 



V[n(a)} 



N r n(a) 
N s {2y^a) d ' 



After some calculation, this result in 

Thus, the operation count and memory needed for SDM 
scales like 



operation ~ N s 



memory ~ N s 



1 



(d+4)/2 



(d+4)/2 



(32) 



We can compare the computational efficiency of SDM 
and spectral (bin) method by comparing the exponents 
in (f29|) and (|32|) . The result suggest that the operation 
count of SDM becomes lower than spectral (bin) method 
when the condition 



d> 



4fc 
4- k 



and k < 4, 



is satisfied, and the memory of SDM becomes lower than 
spectral (bin) method when 



d> 



Ak 
2-k 



and k < 2. 



Thus, if we assume that the attribute-spatial accuracy 
of spectral (bin) method k — 1 ~ 2, SDM becomes more 
efficient in operation count when d > 2 ~ 4, and becomes 
efficient in memory when d > 4 ~ oo. This result may be 
improved if we choose a different kernel function W(a) 
for the density estimator h(a). 

Though our discussions in this section are very brief 
and the assertions must be still verified, our estimation 
suggests that the computational cost of SDM becomes 
lower than spectral (bin) method when the number of at- 
tributes d becomes larger than some critical value, which 
may be 2 ~ 4. 



VII. SUMMARY AND CONCLUDING 
REMARKS 

In this paper, we introduced a simple SDM for warm 
rain and discussed the theoretical foundations and char- 
acteristics of SDM. 

The principle model was introduced at the beginning, 
which is a very primitive microphysics-macrodynmics 
coupled cloud model. We saw how the traditional meth- 
ods can be derived from the principle model and what 
are the problems involved in these methods. Then, we 
developed a simple SDM for warm rain, which can be re- 
garded as a coarse-grained model of the principle model. 
Several numerical simulation schemes for SDM were also 
proposed. Especially, a novel Monte Carlo scheme for 
the stochastic coalescence process was developed and val- 
idated by comparing the numerical results with the solu- 
tions of SCE. We demonstrated the practicality of SDM 
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by presenting a very preliminary result of a shallow mar- 
itime cumulus formation initiated by a warm bubble. 
Brief discussions were carried out on the possibility to 
extend SDM to incorporate various properties of clouds, 
such as, several types of ice crystals, several sorts of sol- 
uble/insoluble CCNs, their chemical reactions, electrifi- 
cation, and the breakup of droplets. The assertions must 
be still verified, but our theoretical estimation suggested 
that SDM becomes computationally more efficient than 
spectral (bin) method when the number of attributes d 
becomes larger than some critical value, which may be 
2-4. 

Though several extensions and validations are still nec- 
essary, we expect that SDM would be more feasible for 
the modeling of complicated cloud microphysics, and pro- 
vide us a new approach to the cloud-related open prob- 
lems, such as the cloud and aerosol interactions, the 
cloud-related radiative processes, and the mechanism of 
thunderstorms and lightning. 
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APPENDIX: HOW TO MAKE A RANDOM 
PERMUTATION 

Let ci^™* 1 = (a^ , , . . . , offl ) be a permutation of 
I n := (1,2,..., n). Let S n := {oW} be the set of all the 
permutation of I n . Thus, S n has n\ number of elements. 
Let us define the random permutation of I n by such a 
permutation that all the € S n occurs with the 
same probability, i.e., P(a (n) ) = 1/n!, Va (n) e S n . 

We have two propositions. 1) = (1) is a ran- 
dom permutation of I\. 2) We can make a random per- 
mutation of I n +\ from a random permutation of /„ by 
the following procedure: Let be a random permuta- 
tion of I n and choose m from 1, 2, . . . , n + 1 randomly. 
Let a(™ +1 ) = (oW,n+l) = (<#°, . . . , at\ n + 1). If 
m =f= n + 1, exchange the value at m and n + 1, i.e., 
ai+i 1 " 1 = dm\ am +1 ^ = n + 1. Then, a^ n+1 ^ is a random 
permutation of I n +i- 

Proposition 1) is obvious and 2) is also easy to 
prove. For all a^™ +1 ' e S n +\, we can uniquely de- 
termine a*-™-* G S n which have the possibility to 
become a^ n+1 \ Indeed, such a*-™-* is given by = 

f (n+1) (n+1) (n+1) (n+1) (n+1) (n+1)-, 

\ a l 7 a 2 i---i a m~l ' a n+l ' a m+l a n ), 

here m satisfies at +1 ' = n + 1. If m = n+1, 
af 1 ^ = a|™ +1 \ i = l,...,n. Then, the 
probability that a^ n+1 ^ occurs is evaluated by 
P(o(" +1 )) = P(oW) x l/(n+l) = l/(n+l)!. 

Consequently, we can make a random permutation of 
I n with 0(n) operations. 
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FIG. 3: Shallow maritime cumulus formation initiated by a warm bubble simulated by SDM. This is a snapshot at the time 
1608 s. All the super-droplets are plotted with the color and the alpha transparency which are determined by their radius Ri and 
multiplicity £j. The color map is indicated in the figure and the alpha value is determined by at = 7.8 x 10~ 5 (-Ri/10 -3 m) 2 £i, 
which is proportional to i? 2 £i. We can see that there is a turbulent like structures inside the cloud. 



